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Abstract 



We consider the problem of pricing an American contingent claim whose payoff depends on 
several sources of uncertainty. Using classical assumptions from the Arbitrage Pricing Theory, 
the theoretical price can be computed as the maximum over all possible early exercise strategies 
of the discounted expected cash flows under the modified risk-neutral information process. 

Several efficient numerical techniques exist for pricing American securities depending on one 
or few (up to 3) risk sources. They are either lattice-based techniques or finite difference 
approximations of the Black-Scholes diffusion equation. However, these methods cannot be 
used for high-dimensional problems, since their memory requirement is exponential in the 
number of risk sources. 

In this paper, we present an efficient numerical technique that combines Monte Carlo simulation 
with a particular partitioning method of the underlying assets space, which we call Stratified 
State Aggregation (SSA). Using this technique we can compute accurate approximations of 
prices of American securities with an arbitrary number of underlying assets. Our numerical 
experiments show that the method is practical for pricing American claims depending on up to 
400 risk sources. On all problems for which we could compare the method with known optimal 
solutions, the price computed through stratified state aggregation was indistinguishable from 
the optimal theoretical price. Several numerical examples are presented and discussed. 
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1 Introduction 

Since the seminal work of Black and Scholes (1973) and Merton (1973) in the early 1970s, 
the arbitrage principle underlying option valuation theory has been extended to a broad range 
of other financial instruments (see e.g. Ross (1976), Cox and Rubinstein (1985)). Indeed, 
any security whose returns are contractually related to the returns on some other security or 
group of securities can theoretically be valuated using the same arbitrage principle. In some 
cases, explicit closed form analytical formulas for the computation of the arbitrage price can be 
derived from this theory. In particular, the original paper of Black and Scholes (1973) provides a 
closed form solution for a European option on a single common stock. Unfortunately, few other 
cases can be solved analytically, and computing the arbitrage price often requires numerical 
simulations. Following an idea initially presented in an early edition of Sharpe (1985), Cox 
et al. (1979) developed a discrete model for the valuation of an American option on a single 
stock that can be easily computed numerically. However, the effective implementation of the 
arbitrage principle is not always such an easy task, and may sometimes become intractable. 
Tractable algorithms have been developed recently for pricing European contingent claims 
with many underlying assets (see e.g. Barraquand (1993)). However, these algorithms cannot 
be used for pricing American contingent claims. 

There are several reasons motivating the development of efficient methods for multidimensional 
contingent claim pricing. In particular, applications exist in the pricing of Over The Counter 
(OTC) warrants, path dependent instruments (Barraquand and Pudet (1994)), multidimensional 
interest rate term structure contingent claims (Heath et al. (1992)) such as mortgage-backed 
securities, and life insurance policies (Fabozzi and Pollack (1987)), futures contracts with 
quality delivery options (Cheng (1987); Boyle (1989)). Also, pricing models taking into 
account the stochastic nature of volatility (Wiggins (1987); Dothan (1987); Hull and White 

(1988) ) require multidimensional modeling. Other applications exist in assets and liabilities 
management, and in corporate capital budgeting (see e.g. Mason and Merton (1985); Coppeland 

(1989) ; Brealey and Myers (1991)). Finally, applications exist in property/liability insurance 
(Merton (1977); Smith (1979); Kraus and Ross (1982); Doherty and Garven (1986); Cummins 
(1988); Shimko (1992)). The tremendous development of financial engineering during the past 
decade can be expected to continue, and new types of securities requiring multidimensional 
modeling are likely to appear at a sustained pace in the future. 

For pricing purposes, financial assets can be divided into two majors classes. The first class is 
that of assets whose future cash-flows cannot be influenced by decisions from the holder taken 
after the purchase date. We will call European instruments all financial assets belonging to 
this first class. In particular, stocks, bonds, futures contracts, European options, swaps, caps, 
floors, mortgage-backed securities are European instruments. The second class is that of assets 
whose cash flows can be influenced a posteriori by the holder. American options belong to this 
class. As another example, the crediting policy of a Life Insurance company selling SPDAs 
greatly influences the present value of the liabilities of the company. This crediting policy can 
be adjusted by the company after signature of the contracts with the policyholders. Therefore, 
the liability associated with the sale of SPDA can be viewed as an American security. We will 
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call American instruments all financial assets belonging to this second class. 

Following the general theory of arbitrage pricing, the theoretical price of a European contingent 
claim is the discounted expected value of its future cash flows under the so-called "risk- 
neutral" probability distribution of the underlying economic factors (Harrison and Kreps (1979); 
Harrison and Pliska (1981); Duffle (1988); Karatzas and Shreve (1988)). Mathematically, 
computing the arbitrage price reduces to computing an integral (sum) over the space of the 
underlying economic factors. When the dimension of the space of the underlying economic 
factors is small, standard techniques for numerical integration can be used. In some cases, 
the integral can even be computed analytically (e.g. Black-Scholes formula). However, the 
computational complexity of evaluating the integral is clearly exponential in the dimension of 
the space. Efficient numerical techniques for pricing high-dimensional European claims are 
presented in Barraquand (1993). 

The price of an American claim is the maximum over all possible cash flow monitoring strate- 
gies of the associated present values of cash flows. For example, the price of an American 
option is the maximum over all possible early exercise strategies of the corresponding present 
values. Since the space of cash flow monitoring strategies is generally huge, direct maximiza- 
tion of the present value is rarely practical (see Bossaerts (1989) for a discussion). However, 
when the underlying economy is modeled as a Markov process, one can use the Bellman 
principle of dynamic programming (Bellman (1957)) to compute the optimal monitoring strat- 
egy. American options are typically priced using a discrete approximation of the dynamic 
programming principle. This is the case in particular of the CRR model (Cox et al. (1979)) 
for American stock option pricing. This approach becomes however impractical when the 
underlying economic space has many dimensions, since the dynamic programming algorithm 
requires a memory space exponential in the number of dimensions. This fact is known as the 
"curse of dimensionality" problem for dynamic programming. 

In this paper, we present a particular state space partitioning technique that attempts to circum- 
vent the curse of dimensionality problem for American security pricing. More precisely, we 
partition the space of underlying assets (the state space) into a tractable number of cells, and 
we compute an approximate early exercise strategy that is constant over those cells. The hope 
is that, if the partition is appropriately chosen, the approximate strategy will be close to the 
actual optimal strategy. Such a partitioning technique is classically called a state aggregation 
technique. 

Among the many possible ways of choosing a partition, one solution is to fix a particular 
real-valued function mapping the state (i.e. the prices of the underlying assets) that particularly 
influences the optimal strategy in the problem at hand. We call this function stratification map. 
Then, the partition chosen is a stratification of the state space into thin layers along this map. 
In other words, we limit our search to strategies that only depend upon the stratification map, 
and not upon the entire state itself. We call this technique Stratified State Aggregation. 

In the case of American security pricing, an obvious candidate for the stratification map is the 
payoff of the security, i.e. the function representing the future cash-flows associated with the 
security. When the stratification map chosen is the payoff of the American security, we call 
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the technique Stratified State Aggregation along the Payoff (SSAP). 

After quantization of the payoff, the SSAP method can be combined with Monte Carlo sim- 
ulation techniques in order to compute the set of conditional probabilities corresponding to 
changes in the payoff value over time. Using these conditional probabilities, an approximation 
of the American price can then be computed backwards in time using techniques reminiscent 
from the classical CRR integration method. 

We implemented the SSAP method on American option pricing problems in dimensions ranging 
from 1 to 400. On all problems for which we could compare the SSAP method with known 
optimal solutions, the SSAP price was indistinguishable from the optimal theoretical price. In 
particular, in dimensions 1, 2, and 3, both put and call prices of options on the maximum of 
the underlying assets were computed accurately by the SSAP method. Also, the SSAP price of 
an American call on the maximum of n assets paying no dividends was indistinguishable from 
the European price for n ranging from 1 to 400, in accordance with a well known theoretical 
result 1 . In other cases, no other method exists to compare to our results. However, the SSAP 
price seems to constitute an accurate approximation of the American price in arbitrarily high 
dimensions. 

To the best of our knowledge, the SSAP method is the first capable of computing American 
prices and exercise strategies in high dimensional cases. 

In order to speed-up the Monte Carlo simulation of conditional probabilities, we developed 
an original variance reduction technique called Quadratic Resampling. Quadratic Resampling 
was originally presented in Barraquand (1993) for European security pricing. In this paper, we 
present an extension of the original QR method that applies to both European and American 
asset pricing problems. Quadratic Resampling consists in correcting the samples obtained 
through classical Monte Carlo simulation in such a way that the expected value of any polyno- 
mial of degree two or less in the space variables is computed exactly. Our experiments show 
that QR is very efficient for American pricing problems in up to 10 dimensions. The average 
speedup obtained through QR ranges from 5 to 35, with an average of about 10. In higher 
dimensions (1 1 and higher), the speedup is only of 2 to 3 in average, and slowly decreases with 
the dimension. 

We implemented a parallel version of the SSAP method on a network of workstations equipped 
with a high-bandwidth interconnect (called a workstation farm). We observed a speedup linear 
in the number of workstations in the network. Both measured and simulated parallelization 
experiments are reported in this paper. 

This paper is organized as follows. In Section 2, we relate our contribution to previous work 
in American security pricing, multidimensional asset pricing, and Monte Carlo valuation. In 
Section 3, we recall the usual assumptions on the stochastic processes governing the evolution 
of securities prices, and the main results of the Arbitrage Pricing Theory. In Section 4, we 

1 Indeed, since the discounted prices of assets paying no dividends are martingales under the risk neutral measure, 
the discounted maximum of such n assets prices is a submartingale. Hence, the corresponding European call price 
is always higher than the immediate payoff. 
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briefly review the current numerical methods used for American security pricing. In Section 
5, we present the method of Stratified State Aggregation. In Section 6, we show how SSA 
prices can be computed through Monte Carlo simulation. In Section 7, we present the method 
of Quadratic Resampling. In Section 8, we present numerical experiments illustrating the 
efficiency of the SSAP method. 

2 Relation to other work 

The theoretical analysis of optimal stopping times for early exercise of American options dates 
back to the work of McKean (1965). This theory has then been further developed by several 
authors (Merton (1973); Harrison and Kreps (1979); Bensoussan (1984); Karatzas (1988); 
Jaillet et al. (1988)). Myneni (1992) surveys the theory of American option pricing. 

The most widely used valuation technique for American options with one underlying asset 
is the binomial lattice approach of Cox et al. (1979). Cox and Rubinstein (1985) outline 
the principle of the multidimensional extension. Other numerical valuation techniques are 
presented in Geske and Johnson (1984), Barone-Adesi and Whaley (1987), Barone-Adesi and 
Elliott (1991). 

The valuation of options depending of several underlying assets has been extensively studied. 
Brennan and Schwartz (1979) addresses the problem of pricing options under two sources 
of risk by direct finite-difference approximation of the generalized Black-Scholes equation. 
In this example the two sources of risk are the short term and the long term interest rates. 
The approach is clearly limited to a few assets, since the memory space requirements and the 
computation time are both exponential in the number of underlying assets. Boyle et al. (1989) 
developed a multinomial lattice method for pricing multidimensional options, in the spirit of 
the approach outlined in Cox and Rubinstein (1985). According to the authors, the computation 
becomes very burdensome for more than two assets. In fact, the multinomial lattice approach 
can be viewed as a finite-difference approximation of the generalized Black-Scholes equation 
using an explicit Euler scheme and an appropriate change of variables (Brennan and Schwartz 
(1978)). 

Stulz (1982) presents an analytical solution to the problem of pricing a European option on 
the maximum or minimum of two underlying assets. The analytical solution is generalized in 
Johnson (1987) to the case of an arbitrary number of assets, taking as given the cumulative 
multivariate normal distribution function. Boyle (1989) and Boyle and Tse (1990) developed 
an approximate method for the same problem. Although the problem is solved analytically 
in Johnson (1987), the approximate method does not require preliminary computation of the 
cumulative multivariate normal distribution function. To the best of our knowledge, no closed 
form solutions have ever been obtained for American pricing problems. 

Good reviews of the Monte Carlo method and different variance reduction techniques such as 
antithetic variables, covariates, stratified sampling, importance sampling can be found in many 
sources such as Hammersley and Hanscomb (1964), Zaremba (1968), Haber (1970), Kalos and 
Whitlock (1986), and references thereof. 
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The application of the Monte Carlo method to option pricing was first presented in Boyle 
(1977), in the context of claims contingent to a single underlying asset. It has then been 
used by several authors for the valuation of path dependent contingent claims. In particular, 
the method has been used for pricing mortgage-backed securities (see Schwartz and Torous 
(1989), Hutchinson and Zenios (1991)). Barraquand (1993) presents the method of Quadratic 
Resampling for Monte Carlo valuation of European securities with many underlying assets. 
The Quadratic Resampling method presented in this paper is an extension of this earlier work 
to the American pricing problem. 

3 Arbitrage pricing theory 

The arbitrage pricing theory is described in many textbooks. We refer the reader to Duffie (1992) 
for a comprehensive and rigorous presentation. In the following, we briefly present through 
intuitive arguments the main results of the theory. These developments do not constitute 
mathematical proofs, and are only aimed at illustrating the main concepts underlying the 
computational approaches to arbitrage pricing. 

3.1 Diffusion model of information process 

We model the economy as a finite-dimensional vector of real-valued state variables X(t) = 
(jci (*),..., x n (t) ) , called factors, representing all the information available to investors at time 
t. Since X(t) represents all information available to agents at time t, in a frictionless market, 
prices of securities must be deterministic functions of time and X{t). It is said that securities 
are contingent claims on the state variable X{t). 

For the sake of simplicity, we will assume that the information process X{t) is a diffusion 
process. However, our results on security valuation described in the next sections apply to 
more general types of stochastic processes. If X{t) is a diffusion process, it is a solution of a 
stochastic differential equation of the type (Ito and McKean (1965)): 



The vector A is called the drift of process X. A is the derivative of the expected value of X. 
The matrix T = BB T is the derivative of the covariance of X. The vector W = {w\ , . . . , w„) is 
a n-dimensional standard Brownian motion. 

Often, the variables x,- are prices of securities available on the market, and are therefore strictly 
positive processes. The expected increments and covariance of increments are then expressed 
in relative values: 



dX = A(X(t) 1 t)dt + B(X(t) 1 t 



)dW 



(3.1) 



dxi 



n 



V/G [l,n] 



Xi 



I^Xi {Xl i ■ ■ ■ jX n , t)dt+ ^V,y(xi, . . . 



x ni t)dwj 



with 



H Xi = cii/xi 



by/xt 
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The matrix V = (v(/)(/,/)e[i,n] 2 * s called volatility matrix, and the covariance of relative returns 
is the matrix JC = (kij) {iJ)e[hn] 2: 



JC = W T 



3.2 European securities 



A security is called European security iff future cash flows cannot be influenced by decisions 
from the holder taken after the purchase date (besides of course selling back the security). Then 
the cash flows are only functions of time and information. 

The cash flow generated by a European security C during the time interval dt, assuming C is 
held indefinitely, will be denoted by fc (X(t) , t)dt. 

If the price C(t) of a European security C is positive, we can define the instantaneous relative 
cash flow rate or dividend yield as: 

d c {X{t),t)=fc{X{t),t)/C{t) 

Let nc denote the expected capital rate of return of C : 

Hcdt = E t {—) 

In the above formula, E t denotes the expectation conditional to the information available at 
time t. The expected total rate of return of C, i.e. the expected rate of return of the total gain 
process is: 

Mg c = +d c 

The dividend yield of the money market account £, called short term interest rate, is denoted 
by r(X). If the proceeds of the money market account are continuously reinvested, the total 
gain process L(to, t) follows the equation: 

L(to, t 0 ) = 1, dL = r(X)L(t Q , t)dt 

or equivalently: 

L(t 0 ,t) = exp^ r(X(r))dT^ 



3.3 Arbitrage pricing 



For any risk factor x,-, let e(C/xi) denote the elasticity of price C to x,-: 

(nl x XidC 

e{c/xi) = 

The major result of the Arbitrage Pricing Theory is the following. There exist numbers 
A i , . . . , X n , called market prices for risk, such that for any security C , the following relationship 
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holds: 

n 

HG c = r + ^e{C/xi)\i 

1=1 

Furthermore, if factor x, is traded, its market price for risk is the expected total rate of return 
on Xi in excess of the riskless interest rate. In particular, if allx,- are traded, we have: 

n n 

1=1 1=1 

Using the properties of diffusion processes, the above results lead to the following partial 
differential equation, called Black-Scholes equation: 

dC (J ^dC, A , 1^ d 2 C 

- w-<*- r)c+ £ - A <>* + 2 E <3 - 2) 

(=1 ij 1 J 

If factor x,- is the price of a traded security, the term fi x . - A, is simply r — d Xi . In particular, if all 
factors are traded, the above equation does not depend on the market prices for risk. Therefore, 
we can replace the information process X by the so-called risk-neutral information process for 
which all market prices for risk are zero: 

dx- " 
Vi G [1 , n] , — = (r - 4,.)^ + Vijdwj (3.3) 

Using a theorem known under the name of Feynman-Kac Formula, we can represent explicitly 
the solution of the above equation: 

where E t represents the expectation under the fictitious risk-neutral information process X 
following the equation (3.3). 

3.4 American securities 

A security is called American security iff it is not European, i.e. if future cash flows can be 
influenced by decisions from the holder taken after the purchase date. Then, the cash flows 
are not only functions of time and information, but also functions of the decisions taken by the 
security's holder. A cash flow monitoring strategy u is a stochastic process associating with 
each time and information state a decision u(X(t), t) e U, U being an appropriate decision 
space. We denote by CMS the space of cash flows monitoring strategies. For the sake of 
simplicity, we will consider only finite American securities, i.e. those for which CMS is a finite 
set. The cash-flow generated by an American security C during the time dt can be written: 

fc(u,X, t)dt 
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Let C be an American security with a cash flow function fc (w, X, t). To each cash flow 
monitoring strategy u we can associate the European security C" whose cash flow function/^ 
is: 

f?(X, t)=f c (u(X, t),X,t) 

Let us assume that the market is complete, i.e. that any factor x,- can be perfectly hedged by 
building an appropriate dynamic portfolio of traded assets, called trading strategy. Then, for 
any cash flow monitoring strategy u, the security C u can be replicated by a trading strategy. 
We can therefore consider its arbitrage price C u . In particular, we can consider the European 
security C u * maximizing the market value: 



Barring arbitrage, we must have: 



C" = max C" 
»eCMS 



C M * = max C u 



«eCMS 

Indeed, if C < C u * , we can buy C, sell C u * , and take the proceeds C u * — C immediately. Then, 
by selecting the cash flow monitoring strategy u* for C, all future cash outflows generated by 
the sale of C u * will be exactly compensated by the cash inflows generated by C, so that we 
will not be obligated to any future payments. In effect, we will have made money without any 
investment and without risk. 

Reciproqually, if C > C u * , let u$ be the cash flow monitoring strategy chosen by a purchaser 
of security C. Since by definition of u* we have C u > C u ", we must have C > C"°. By selling 
C and buying C u ", the buyer can immediately take the proceeds C - C"°, without any change 
in the future cash flows. This is again an arbitrage opportunity. 

We can state: 

In a complete market, the price of an American security is the maximum over all 
possible cashflow monitoring strategies of the corresponding European prices. 

In other words, computing the price of an American security C reduces to computing that of 
the equivalent European security C u * . 

C(X( t) , ,) = max E, ( f f^hlh^Il dT 



«eCMS V it L \f T ) 

Therefore, the differentiation between European and American securities is irrelevant in fi- 
nancial terms in a complete market free of arbitrage opportunities. On the other hand, if the 
market is not complete, there is no specific relationship between the American price and the 
corresponding European prices. 

As a computational problem however, American security pricing is much harder that European 
security pricing. Indeed, the determination of the equivalent European security requires to 
precompute the optimal cash flow monitoring strategy. In many cases, this precomputation 
step is practically much more complex than the European price computation itself. 
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3.5 European and American options 
3.5. 1 European options 

We consider an arbitrary asset S that can be replicated by an appropriate trading strategy. In a 
complete market, we can assume that the information vector X = (jci , . . . , x n ) represents the 
prices of n given traded securities. 

For example, S could be the right to purchase any one of the n securities at a given expiration 
date T. Barring arbitrage, the price S(X(T), T) at expiration date is 

S(X(T),T) = max jc,-(T) 

ie[l,n] 

More generally, the price 5 can be any contractual function g of the underlying securities prices 
x\ , . . . , x n at expiration. 

S(X(T),T) = g(x 1 (T),...,x n (T)) 

The function g, which represents the unique cash-flow associated with the contractual asset S, 
is called the payoff function. 

By definition, a European call (resp. put) option on an asset S with expiration date T and strike 
price K gives its owner the right to purchase (resp. sell) at time T the asset S for the price K. 
Since they leave a choice to the owner, European options are theoretically American securities. 
Indeed, the holder can choose to exercise or not exercise the option at expiration date. The 
decision space corresponds to this choice: U = {exercise, no-exercise}. 

The cash flow monitoring strategy is any process associating the decision to the time t and the 
available information at t. The space of admissible cash flow monitoring strategies CMS is 
composed of all adapted processes u verifying: 

yt<T,yX, h(X, t) = no-exercise 

and taking the two possible values at time T: 

u(X, T) G {exercise, no-exercise} 

For the sake of simplicity, we will assume that the payoff of the option at exercise time T is 
distributed during a small time interval [T, T + At] . We will also assume that exercise decisions 
u(X(t) , t) are piecewise-constant, i.e. change only at the beginning of time intervals of duration 
At. For notational convenience, we will define the "spike" function 8{u) associating the value 
1/A? to the decision exercise and the value 0 to the decision no-exercise. 

The cash flow function of a call option is: 

f{u,X,t) = {S{X,t)-K)5{u) 
The optimal exercise strategy u* is easily identified: 

Mt^T, ?) = no-exercise 
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and 



ii* (X, T) 



exercise if S(X, T) > K 



no-exercise otherwise 



Therefore, we identify the call option with its canonically associated optimal European security, 
whose cash flow function is: 



Hence, from a computational viewpoint, the European call option can be considered as a 
European security 

If the payoff at expiration date T is g(X), the price of the European call option can be written: 



More generally, any European contingent claim with a single cash flow / at date T can be 
priced according to the above formula. Efficient numerical techniques based upon Monte 
Carlo simulation exist for computing numerically the above expectation in arbitrarily many 
dimensions (see e.g. Barraquand (1993)). 

3.5.2 American options 

By definition, an American call (resp. put) option on an asset S with expiration date T and 
strike price K gives its owner the right to purchase (resp. sell) on or before time T the asset S 
for the price K. The space of admissible cash flow monitoring strategies CMS is composed of 
all adapted processes u taking the two possible values: 



and verifying: 

Mt < T, u(X(t), t) = exercise => Vr > t, u(X(t), t) = no-exercise 

That is, exercise cannot occur twice. The cash flow function of an American call option is 
identical to that of a European call option. Unfortunately, in this case, the computation of 
the optimal early exercise strategy is not straightforward, since exercise can happen before 
expiration. We can simplify the above formulation by noticing that if it is optimal to exercise at 
a given underlying asset price So, it is also optimal to exercise at any higher price. Therefore, 
if we denote H(X(t), t) the smallest possible value of So, the optimal early exercise stopping 
time r* is the solution of the following equation: 



f*(X, t) =f(u*(X, t),X, t) = max(S(X, t) - K, 0)S(u*) 




with 



f(X) =max(0 1 g(X)-K) 



Mt G [0, T], u(X, t) G {exercise, no-exercise} 



t* = inf {/ G R+ , S(X(t) , t) = H(X(t) , t 



")} 
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The following arbitrage argument shows that 

H = K+C 

where C is the American call option price. Indeed, whenever C < S - K, any investor could 
buy the option, exercise it immediately, and take the proceeds (S - K) - C > 0. Therefore, it 
is optimal to exercise whenever 5 > K + C, hence H < K + C. 

Reciprocally, if C > S - K, no investor holding the option would be willing to exercise it and 
take S - K, since by just selling it he would make a higher immediate profit C > S - K. At 
any time, and for any information state, the optimal stopping time r* verifies: 

C(X(r*), t*) = S{X(t*),t*) - K 

Hence, the optimal early exercise strategy can be written: 

u *(Xt)-i exercise if C(X,t) <S(X,t)-K 
' \ no-exercise otherwise 

We see that the computation of the optimal early exercise strategy requires to precompute the 
pricing relationship between the option and the underlying asset, which is what we were trying 
to compute in the first place. 

If the payoff at the exercise date is g(X), the price of the American call (resp. put) option can 
be written: 

C(X(t),t)= max £ f (^ffllll5( M (X(r),r))A/) (3.4) 

i.(x(T),T),T>(,«eCMS L{t,T) J 

with 

f(X) = maL(0,g(X)-K) (resp. max(0, K - g(X))) 

More generally, any contingent claim entitling its holder to the single cash flow/ on or before 
an expiration date T can be priced according to the above formula. 



4 Numerical methods for American security pricing 
4.1 Stochastic dynamic programming 

The explicit numerical valuation of an American option using the above formula involves a 
maximization over the set of of all possible early exercise strategies. The strategy u can be 
any function associating to each current value of the underlying assets X = (xi, . . . ,x n ) and 
each current time t the binary decision to exercise or not exercise. Since the number of such 
possible strategies is huge, direct maximization is rarely practical (see Bossaerts (1989) for a 
discussion). 

The only practical technique to date consists in using the Bellman principle of Dynamic 
Programming. This principle can be applied since the information vector X is assumed to be a 
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Markov process, and therefore the optimal early exercise strategy u(X(t) , t) only depends upon 
time and the current vector X(t). 

Assuming that exercise decisions can only be taken at discrete times intervals of constant dura- 
tion At, the maximization problem (3.4) can be rewritten using the law of iterated expectations 
and the property L(t, f ) = L(t, t + At)L(t + At, f ) : 



C{X(t),t)= max E, 
n(x(f),f),weCMS 



f{X ^ 8(u(X(t),t))At+ 



L(t, t + At) 

T 



max Et+Ati uT+Kt J ^ X{T ^ T ^ 

u(X(T),T),T>t+At,ueCMS Vr^+Af ^ + ' > J 

Examining successively the two cases u(X(t), t) = exercise and u(X(t), t) = 
no-exercise, and using the expression of C(X(t + At), t + At) from equation (3.4), we 
obtain 

C(X(t),t) = exp(-rAt) max (f(X(t)),E t (C(X(t + At),t+ At))) (4.1) 

assuming that the interest rate r is piecewise constant on intervals of duration At. The above 
recursive expression, called Bellman equation, allows to compute the price C of an American 
option by proceedings backwards in time from the expiration date T. Using the properties 
of diffusion processes (Ito's formula), it can be shown (see e.g. Jaillet et al. (1988)) that 
the solution of the above equation (4.1) converges towards the solution of the Black-Scholes 
equation (3.2) when At converges towards 0. 

dC ^dC, , , lr. d 2 C , 

" -Yt = ~ rC + 2. ^ r - + 2 2. 

(=1 ij 3 

with the additional time-dependent boundary condition: 

C(X,t) >f(X) 



4.2 Finite differences and the Cox-Ross-Rubinstein approach 

The general method for solving the above partial differential equation (PDE) is to quantize it 
using a finite difference method (see e.g. Duffle (1992), chapter 10). 

We will illustrate this approach using a number of simplifying assumptions. We will assume that 
the process of the underlying securities is jointly lognormal, i.e. that the mean and covariance 
matrix of relative returns are constant. We will also assume that the interest rate r and the 
dividend yields of the underlying securities are constant. Then, the change of variable below 
simplifies the finite difference approximation. We consider the vector Y = (yi , . . . , y n ) T - 

Vi G [1, n], y t = logXi - (r - d Xi - ^k ti )t 

and we then define the vector W = (w\ , . . . , w n ) T 

W= V~ l Y 
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By construction, Y = VW, and W follows a k-dimensional standard Brownian motion. The 
Black-Scholes equation in the variable W writes: 

_ rt d{e- rt C) _lsr- d 2 C 

~ 6 Ft ~2^dwj 
i 1 

Defining 

= (wi, . . . , w,- + Aw,-, . . . , w n ) T , = (wi, . . . , Wi - Aw,-, . . . , w n ) T 

then writing 



<PC _ C{Wf,l + Al) - 2CiW r + At) + CjWr , , + AQ ^ 
awf Awf 

and 



V(e-*C) _ g ^C(W,Q-C(W,f + Af) 

^ " A/ + 1 j 

with 

Vi G [1, n], Aw,- = a/«a7 

we get the simple explicit Euler scheme (see Duffie (1992) for a description of more sophisti- 
cated schemes): 

1 " 

C(W, t) = e~ rAt - C{Wt ,t + At) + C(Wr , t + At) 

i=l 

with the terminal boundary condition: 

C(W(T) 1 T)=f(X(T)) 

and the early exercise condition: 

C(W(t),t)>f(X(t)) 
where X(t) = (xi(t) , . . . , x„(f)) is obtained by the inverse formula: 



/ 1 n 

Vi G [1, n], x,-(f) = x,-(0) exp (r - 4 ; - -k u )t+ ^ v ;y w fc (?) 

V 



(4.2) 



We implemented the above numerical scheme for an arbitrary number n of underlying assets. 
The method yields accurate results for n < 3, but its memory requirement is intrinsically 
exponential in n. Hence, it cannot be used for n > 3 with a quantization step At small enough 
to yield accurate results (See section 8). Many variants and generalizations of the above 
finite-difference method have been studied by a number of authors (see Duffie (1992)). 

An alternative technique consists in quantizing directly the Bellman equation (4.1). This is 
the paradigm underlying the original Cox-Ross-Rubinstein approach (Cox et al. (1979)). The 



Research Report No. 38 



April 1994 



14 



Jerome Barraquand and Didier Martineau 



Brownian motion W is approximated by an w-dimensional binomial process W At defined as 
follows: 

Ve= (ei,...,e„) G {-1,1}", Prob(W At (t + At) = W At (t) + eVAt) = 1 
Then, the quantized Bellman equation writes: 

C(W,t) = exp(-rAt)m a xlf(X(t)),^ £ C(W + eVAt, t + At)) (4.3) 

\ =€{-1,1}- / 

For n = 1 , the two formulations are exactly equivalent. In higher dimensions, they yield almost 
identical results for small enough At. In summary, the finite-difference and binomial lattice 
methods are essentially equivalent, and both intractable for n > 3 due to their exponential 
memory requirement. 



5 State aggregation 

5.1 State aggregation price 

Classical numerical methods being unable to deal with high-dimensional American valuation 
problems, one must resort to alternate approximation schemes. State aggregation is a classical 
approximation technique for the numerical solution of stochastic optimal control problems (see 
e.g. Bertsekas (1987), Kushner and Dupuis (1992)). 

For the problem of American security pricing, the relevant state space in the Bellman equation 
(4.1) is the ^-dimensional space of the underlying assets values X = (xi, . . .,*„). State 
aggregation consists in partitioning the state space into a tractable number of cells, and in 
computing an approximate early exercise strategy u(X, t) that is constant over those cells. The 
hope is that, if the partition is appropriately chosen, the approximate strategy will be close to 
the actual optimal strategy. 

Mt G {0, At, . . ., T}, let us consider a finite partition P(t) = (P\(t), . . .,Pfc(f)W) °f tne state 
space R n + , i.e. a set of k(t) subsets of R n + verifying: 

|J Pi{t) = R\ and V(U) G [l,k(t)] 2 , i^j, P^t) f)Pj(t) = 0 

se[i,*(0] 

We assume that the partition P(0) only has two cells: 

Pi(0) = {X(0)}, and P 2 (0) = R n + \{X(0)} 

Among the set CMS of all possible early exercise strategies, we consider the subset U{P) of 
piecewise constant strategies, i.e. of strategies u(X, t) that are constant over each cell Pj(t) of 
the partition. 
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Then, we define the state aggregation price C* (i, t) as the maximum over all possible piecewise 
constant strategies in U (P) of the expected risk-neutral discounted future cash-flow condition- 
ally to the event X(t) G Pi{t) : 



C*(M) = max £ 0 ( E^y^WM, r))Af | *(*) G P,(0 J (5-D 

u(X(r),T),T>t, ueU(P) \^ L(t,T) J 

Since Pi(0) = {X(0)}, and since U{P) C CMS, the state aggregation price at initial time 
C*(1,0) is obviously upper bounded by the true American price C(X(0),0). Furthermore, 
since the strategy ue U w corresponding to the European price consists in never exercising before 
the expiration date, it is clearly constant over the cells of any partition before expiration. Hence, 
by definition of the state aggregation price, the European price is upper bounded by the state 
aggregation price at the initial time C* ( 1 , 0) . We can state 

For any family of finite partitions P, the state aggregation price C* = C* ( 1 , 0) is 
lower bounded by the European price and upper bounded by the American price. 

qEuto ^ q* ^ QAmer 

We will now derive a recursive backward valuation formula for the state aggregation price, 
under an additional Markovian assumption. 

5.2 Markovian approximation 

We will now assume that the partition P is such that the process l(t) defined by X(t) G Pi( t ) if) 
is approximately Markov under the risk neutral measure: 

Vi,j,0, E Q (<f>(X(t+2At)) | X(t) G Pi(t), X(t + At) G Pj(t + At)) 
« E 0 {<f>{X{t + 2At)) | X(t + At) G Pj(t + At)) 

Applying again the law of iterated expectations to the definition of the state aggregation price 
in equation (5.1), we get: 

f{X ^ 8(u(X(t),t))At+ max 



L(t, t + At) u{X(r),T),T>t+At,ueU(P) 



C*(I(t),t)= max E Q 
u(x(t),t),ueu(p) 

Eo \ E Z^+StT ^ m(x(t) ' T))At 1 X{t) e Pl ^ th ^ + Af )eP/( f+ Ao(^ + AO] 

\mep m (t)] 

By examining successively the two cases u(X(t), t) = exercise and u(X(t), t) = no-exercise, 
and applying formula (5.1) at time t + At: 

CT(/(0,0« 

e-^max {E 0 (f(X(t)) \ X(t) G P I(1) (t)), E 0 {C*(I(t + At), t + At) \ X(t) G P m (t)))) 
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where the substitution of formula (5.1) at time t + At is justified by the above Markovian 
approximation. 

By definition of a partition, we have: 

k(t+At) 

1 = Y lx(t+At)EPj(t+At) 

where 1^ is the indicator variable of the event A (i.e. takes the value 1 iff A is true, and 0 
otherwise). Combining this with the above recursive equation, we get: 

/ k(t+At) 

Vie[l,*(f)], C*(i,t)^e- rAt m a x U(t), Pij(t)C*(j,t + At) 
where we have defined for notational convenience 

f i (t)=E 0 (f(x(t))\x(t)ePi(t)) 

and 

Pi j(t) = Prob (X(t + At) e Pj(t + At) | X(t) e Pi(t)) 
Furthermore, the value at expiration date is determined by the terminal condition: 

C*(i, T)=MT) 

5.3 Recursive state aggregation 

We define the recursive state aggregation price Csa as the solution of the following program: 

C SA {h T)=fi(T) 

and 

/ k(t+At) 

C SA (i, t) = e~ rAt max U(t), Pij{t)C SA {j,t + At) 

When the partitions P{t) are chosen in such a way that the process I(t) is actually Markovian, 
the recursive state aggregation price Csa is exactly the true state aggregation price C* . The 
implementation of a recursive state aggregation program proceeds in two steps. 

1) Definition of an appropriate family of partitions, 

2) Computation of the expected payoffs fj(t) and the conditional probability matrices py(t). 

Then, the approximate price of an American contingent claim with terminal payoff / can be 
computed backwards in time using the above recursive equation. 
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Finite difference or lattice-based methods can be viewed as particular instances of the recursive 
state aggregation method. Indeed, each grid point or lattice point W can be viewed as the 
center of the hypercube Hyper(W) defined by: 

Hyper(W') = {Z G R n , V/ G [1, n], w) - l -A Wj < Zj < w) + l -A Wj } 

Although the process I(t) such that W(t) G Hyper(Wj( f )) is not actually Markovian, it can be 
approximated by a Markov process with reasonable accuracy for a fine enough lattice grid. In 
the lattice approach, this corresponds to approximating the Brownian motion W by the binomial 
process W Af . Then, the expected value over the hypercube fi(t) is approximated by the value 
in the center/(X(W)). Similarly, the conditional probability is simply taken from that of the 
binomial approximation, which yields the recursive formula (4.3). As explained in section 4.2, 
the limitation of this approach is the exponential growth of the number of lattice hypercubes in 
the number n of underlying assets. 

One solution for generating a partition P(t) with a tractable number of cells is to fix a priori 
a small number k and set for all t > 0, k(t) = k. Then, one can sample for each time t the 
state space with k samples following the risk-neutral distribution of X(t). For example, if X(t) 
follows a jointly lognormal distribution, the samples X' (t) can be computed from the samples 
W through formula (4.2). In turn, the samples W(t) can be generated following a jointly 
standard normal distribution with variance t along each coordinate. Then, the cell Pi{t) of the 
partition P(t) can be defined, in the Brownian Motion space (i.e. the space of the variable W), 
as the set of points W closer to W'(t) than to any other sample. Such a partition is commonly 
called the Voronoi partition associated with the samples W l (t), . . ., W k (t). Unfortunately, 
Voronoi partitions have an undesirable asymptotic property for large n. The cells tend to 
become so large that the probability Pij{t) of moving into cell j at time t + At from cell i at time 
t is almost 0 for all fs but one. Hence, a Voronoi partition of the ^-dimensional state space 
does not accurately reproduce the diffusion effect in the information vector X(t). In order to 
devise partitions such that several of the conditional probabilities py(t) are non-zero, one must 
build cells having a small directional diameter in the Brownian Motion space, i.e. such that 
the numbers 

D, = sup inf 

XePi(t) S=(si,...,S n ),s\+...,S 2 n = \ 

{A max - A min , A min = inf{A G R,X+ \s G P,-(f)}, A max = sup{A G R,X+ \s G P,-(f)}} 

are as small as possible. Indeed, the smaller D, is, the higher is the probability that the Brownian 
motion process W(t) crosses the boundary of the cell Pi(t). 

5.4 Stratified state aggregation (SSA) 

One solution for choosing a partition with small cell directional diameters is to fix a real-valued 
function mapping the state that particularly influences the optimal strategy in the problem at 
hand. We call such a function a stratification map. Then, the partition chosen is a stratification 
of the state space into thin layers along this map. If the layers are chosen sufficiently thin, the 



Research Report No. 38 



April 1994 



18 



Jerome Barraquand and Didier Martineau 



diameters of the cells along the direction of the gradient of the map will be small. Hence, the 
probability of crossing the boundary of a cell during a small time interval will be relatively 
high, and the drawback of Voronoi-based partitions will be avoided. 

In order words, stratification consists in limiting the search to strategies that only depend upon 
the stratification map, and not upon the entire state itself. We call this technique Stratified State 
Aggregation. 

In general, we can consider a vector- valued /-dimensional stratification map (/ < n): 

h: R n xR — > R l (5.2) 
(X, t) — ► h{X, t) 

and a family Q of partitions of R l . From the family Q and the map h, we can build the reciprocal 
image partition P = h~ l (Q): 

Pi{t) = {XeR n + ,h{x,t)eQi{t)} 

In the case of American security pricing, an obvious candidate for the stratification map is 
the payoff of the security. When the stratification map chosen is the payoff of the American 
security, we call the technique Stratified State Aggregation along the Payoff (SSAP). 

Let us consider an American security with a single cash-flow/ (X) on or before an expiration date 
T. In particular, in the case of an American call option, the cash-flow is/(X) = max(0, g(X) - 
K) , where K is the strike price. In order to illustrate the SSAP method, we set / = 1 and choose 
h(X, t) =f(X) for the stratification map. We also take k(t) = k constant for all times t £ [0, T] . 
In the numerical examples developed in section 8, we assume that g(X) = max,^ „] x;, and 
that the process X is jointly lognormal, of the form described in section 4.2. The partitions 
Q(t) of the image space R l = R are chosen logarithmic in all these examples, i.e. the interval 
Qi{t) , t > 0 is of the form: 

V/ e[2,k- 1], Qi{t) = (A(t)e B ^ i - 2 \A(t)e B ^ i -^ 

and 

Gi(f) = (-oo,A(f)] , Qk(t) = (A(t)e B ^ k - 2 \+^) 
for adequate parameters A (f) and B(t) . The cell Pi (t) is then by definition: 

V/ e[2,k- 1], Pi{t) = {Xe R n + , A(f)eW" 2 ) <f(X) < A(f)e B W(''- 1 )} 

and 

p l (t) = {x e R n + , f(x) < A(t)}, p k (t) = {x e R n + , f(x) > A(^W" 2 )} 

In our experiments, the number k of cells is set to 100. The numbers A(t) and B(t) are 
automatically adjusted so as to ensure: 

Prob (X(t) e Pi(t)) ^ Prob (X(t) e P k (t)) « 0.1% 

The numerical results obtained with the SSAP method, presented in section (8), show that these 
empirical parameters are adequate for a broad range of American security pricing problems. 
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6 Monte Carlo estimation of American price 

6.1 Generation of sample paths 

Once the family of partitions P has been chosen, for example using the SSAP method, it remains 
to compute numerically the expected payoffs fi(t) and conditional probabilities Pij(t). These 
numbers can be expressed as integrals over the state space. In general, they must be computed 
numerically. The only general tractable method for computing such high-dimensional integrals 
is the Monte Carlo method. 

It consists in generating a given number M of sample paths for the underlying assets price 
process X(t). In general, this can be done through direct numerical integration of the Ito 
equation (3.3). A simple explicit Euler scheme is given by: 

Xi(t + At) = Xi (t) exp ^jr - d Xi - ^k^j (X(t) , t)At + ^ v^(X(t) , t)VAt z) 

where zj follow independent standard normal distributions for all j and t. d = T/At being the 
number of time steps in [0, T], we must draw a total of M x d X n standard normal variates in 
order to generate M n-dimensional sample paths X 1 (t) , . . . , X M (t) for all t > 0. 

In general, At must be chosen small enough so as to reach a reasonable accuracy. In practice, 
a number of time steps d = 100 is sufficient in most asset pricing applications. However, 
when the joint process X(t) is assumed lognormal as in section (4.2), d can be chosen much 
smaller. Indeed, the underlying assets price process X can then be obtained by formula (4.2) 
from a standard Brownian motion W. In our experiments, we found that a number of time 
steps d = 10 is sufficient for American security pricing with lognormal underlying assets price 
processes. 

6.2 Conditional probabilities and payoff expectations 

Once the M sample paths X 1 (t) , . . . , X M (t) are computed, the number a,-(f) of samples crossing 
Pi(t) and the number &//(*) of samples moving from Pj{t) to Pj(t + At) are easily computed: 

fli(f) = Caxd{ke[l,M],X k (t) ePi(t)} 

bij(t) = Card{k e [I, M], X k (t) ePi(t) and X k (t + At) £ Pj(t + At)} 
Similarly, the sum Cj(t) over of samples X k of payoff values/ (X k (t)) is computed from: 

Ci (t) = /(**«) 

{ke[i,M],xk(t)ePi(t)} 

By the law of large numbers, we have the following identities: 

/ .\ ,■ bij{t) d(t] 



m — s-cxd aAt) m — s-cxd aAt) 
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6.3 Backward integration algorithm 

Using the above Monte-Carlo estimates of the conditional probabilities and payoff expectations, 
an approximation of the American price can then be computed backwards in time using the 
simple algorithm described below. 

• At time T, the approximate SSAP price is initialized at: 



C(i, T) = 




• At time T - At, we can compute for all i £ [!,£]: 



• The above procedure is then applied recursively, backwards in time, to compute all the 
prices C(i,T - 2At),C(i,T - 3At), . . . , C(l, 0) = C SSAP . 

The memory required in the SSAP method is proportional to k 2 x d, corresponding to the storage 
of the conditional probabilities py(t) . The computation time is proportional toMxn 2 x d+k 2 x d, 
the first term corresponding to the drawing of the M Monte Carlo sample paths, and the second 
to the backwards integration. Hence the memory and time complexities of the SSAP method 
are polynomial in n. This is to be contrasted with classical PDE methods which are exponential 
in n. 

7 Quadratic resampling 

7.1 Quadratic resampling for multidimensional Monte Carlo integration 

Each standard normal variable Zj is simulated by generating M standard normal deviates 
zj(l), • • .,Zj(M). Many variance reduction techniques exist to improve the computational 
efficiency of the Monte Carlo method. In particular, we use the well known technique of 
Antithetic Variables (AV), which consists in generating only M/2 standard normal deviates 
zj ( 1 ) , • • • , Zj(M/2) and take for the remaining samples zj(M/2 + 1), . . . , zj(Af) the opposite 
values — zj(l), . . . , — zj(M/2) (M is assumed even). 

In addition to the technique of Antithetic Variables, we also use an extension of the Quadratic 
Resampling technique. Quadratic Resampling was first presented in Barraquand (1993) for 
reducing the variance of multivariate Monte Carlo integration. 

We briefly present below the original QR method before describing the extension we developed 
for the purpose of the SSA method. We consider the problem of computing the integral: 




E(f(X))= [ f(X)p(X)dX 



Jr n 
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where p is a given probability density function over R N , and/ any integrable function. We 
assume that the expected vector E(X) and the covariance matrix Kx have been precomputed. 
In the application we consider here, p (X) is the density of the ^-dimensional standard normal 
distribution. Hence E(X) = 0 and Kx = In- 



Then, we can approximate the integral E(f(X)) by: 

i M 

i=i 

In particular, the empirical mean is: 

M 



(=1 

We can likewise define the empirical covariance: 



Kx~ = (X - X)(X - X) T = XX T -XT 
We define the gain matrix (for M large enough, Kx is regular): 



H = \/K x y (Kx) 

and the new random variable: 

Y = H(X -X) + E(X) 
We consider the M samples Y l = H(X' - X) + E(X). For this particular sampling, we get: 

M 



Similarly 



Y = - > Y l = E(X) 
M 



Ky = (Y — Y)(Y — Y) T = (H(X - X))(H(X - X)) T = HK X H T 
Using the definition of H we get: 

Ky = Kx 

Hence, the empirical first and second order moments using the samples Y l are exactly equal 
to the real first and second order moments of X. In particular, the empirical mean of any 
polynomial/ of degree two or less in the variables jci, . . . ,Xn verifies: 



i M 

/(r) = ^E/( yi ') =i W)) 



!=i 

The method of quadratic resampling consists in using the samples Y l in place of the samples 
X' in the quadrature formula. We can state: 

Any numerical quadrature formula generated through quadratic resampling is 
exact for any polynomial of degree two or less 
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7.2 Quadratic resampling in space Iime 

In the American pricing problem, we must sample not only the underlying assets prices X(T) 
at expiration date, but also those prices X(t) for all possible early exercise dates. Hence, 
quadratic resampling must be applied to the underlying assets space domain elevated to the 
power of the time domain. In other words, we consider the variable Z = (^)je[i,n],te[At,T] in the 
n x (i-dimensional space time domain R" xd . Z is an x <i-dimensional standard normal variable. 
Given M standard normal deviates zj(l), . . . , zj(M) (with antithetic variables) for each of the 
n x d variables zj, we can consider them as M vector samples Z l , . . . , Z M inR" xd . Then, we 
can apply Quadratic Resampling to the n x d- vector variable Z. 

We have E(Z) = 0 and Kz = hxd- Since the samples are generated using antithetic variables, 
the empirical mean Z is 0. We can compute the empirical (n X d) X (n X d) covariance matrix 
Kz- Then, we can apply to each sample Z l the transform: 



Finally, we can replace in the Monte Carlo simulations presented in section (6) the samples 
Zj(i) by the samples components of the n x d- vector y l . Experimental results reported 
in section (8) demonstrate the efficiency of the extended quadratic resampling method for 
American option pricing. 

8 Experimental results 
8.1 A test case 

We implemented the method of Stratified State Aggregation along the Payoff function (SSAP). 
We also used the quadratic resampling technique for drawing the Monte Carlo sample paths. 
We present below some numerical results for several European and American option pricing 
problems ranging from 1 to 400 underlying assets. In all the experiments, we assumed that 
the American options can only be exercised at d = 10 different dates during the life of the 
option. This corresponds to choosing a time step At = T/d = T/10. We experimented with 
several different payoff functions, in particular with payoffs corresponding to the maximum, 
the minimum, or the average of the n underlying assets. We obtained similar results for all 
these different payoff functions. We only present below the case of an option on the maximum 
of the underlying assets. Its payoff function is defined by: 



where K is the strike price of the option. 

In all the experiments, we assumed that the underlying assets price process X(t) is lognormal, 
with a covariance matrix of relative returns JC of the form: 




g{x\, . . . 



x n ) = max(0, max(xi, . . . 




Vi G [l,n], ku 
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Call option prices: x\ (0) = 


$40, r = 


5% 


Parameters 


European 


American 


c 1 


T 
1 


A 




^ssap "rstaev 




*-SSAP 


Astdcv 






JJ 


J.1 J 


J.1 J U.UU1 


J.lJ 




u.uuz 




1 


40 


1.00 


1 00 0 006 


1.00 


1.00 


0.008 








U.UZ 




u.uz 


n no 

\jXjL 


U.UU1 






JJ 


J. /O 


J./O U.UUJ 


J. /o 


J. /o 


U.UUJ 


20 % 


4 


40 


2.16 


2 16 0 006 


2.16 


2.16 


0.008 






T"J 


UJU 




U.JU 


UJU 


n nnzL 






JJ 




f. a(\ n c\r\A 






u.uuo 




7 


40 


3.00 


2.99 0.010 


3.00 


3.00 


0.010 






45 


1.09 


1.09 0.006 


1.09 


1.09 


0.006 






35 


5.38 


5.38 0.003 


5.38 


5.40 


0.007 




1 


40 


1.91 


1.91 0.010 


1.91 


1.92 


0.010 






45 


0.41 


0.41 0.003 


0.41 


0.41 


0.003 






35 


6.88 


6.88 0.010 


6.88 


6.90 


0.020 


40% 


4 


40 


3.96 


3.96 0.020 


3.96 


3.97 


0.020 






45 


2.08 


2.08 0.006 


2.08 


2.09 


0.010 






35 


8.07 


8.08 0.010 


8.07 


8.10 


0.020 




7 


40 


5.35 


5.34 0.024 


5.35 


5.36 


0.040 






45 


3.40 


3.40 0.012 


3.40 


3.42 


0.020 





Put option prices: . 


d(0) =$40, r = 5% 


Parameters 


European 


American 




T 




Pbs 


Pssap Astdev 


Ppde 


Pssap 


Astdev 






35 


0.00 


0.00 0.001 


0.00 


0.00 


0.001 




1 


40 


0.83 


0.83 0.008 


0.84 


0.84 


0.008 






45 


4.84 


4.84 0.001 


5.00 


5.00 


0.000 






35 


0.19 


0.19 0.003 


0.19 


0.19 


0.004 


20% 


4 


40 


1.50 


1.50 0.006 


1.56 


1.56 


0.010 






45 


4.77 


4.77 0.004 


5.06 


5.07 


0.010 






35 


0.41 


0.41 0.003 


0.42 


0.42 


0.006 




7 


40 


1.86 


1.86 0.010 


1.96 


1.96 


0.012 






45 


4.82 


4.82 0.006 


5.24 


5.23 


0.016 






35 


0.24 


0.24 0.003 


0.24 


0.24 


0.003 




1 


40 


1.74 


1.74 0.010 


1.75 


1.76 


0.010 






45 


5.23 


5.23 0.003 


5.27 


5.29 


0.015 






35 


1.31 


1.31 0.010 


1.32 


1.33 


0.020 


40% 


4 


40 


3.31 


3.31 0.020 


3.36 


3.37 


0.020 






45 


6.35 


6.35 0.006 


6.47 


6.49 


0.015 






35 


2.08 


2.09 0.012 


2.12 


2.13 


0.015 




7 


40 


4.21 


4.21 0.020 


4.31 


4.31 


0.016 






45 


7.12 


7.13 0.012 


7.36 


7.34 


0.030 



Table 1: Results of the SSAP method with 1 underlying asset 
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and 

V(U) G [1, n] 2 , i 7^ J) fcy = /OCTjtfy 
for « + 1 numbers a\ > 0, . . . , a n > 0 and — l/(n — 1) < p < 1. 

Volatilities (cr, ), correlations (p), and interest rate (r) are counted in percent per year. The time 
to expiration T is counted in months, with the convention 1 month = 30 days. All asset and 
strike prices are counted in dollars. 

Since SSAP uses Monte Carlo simulation, we report confidence intervals together with all 
results. These confidence intervals are computed from the central limit theorem, i.e. we 
assume that at a confidence level of 99.95%, the error must be less than 4 times the observed 
standard deviation of the result. The confidence interval reported is 4 x stdev. 

All the simulations were run on a DEC 3000 model 500X workstation, with an ALPHA AXP 
processor running at a clock rate of 200 Mhz, and 1 Gigabyte of main memory. 

8.2 One underlying asset 

We first study the one-dimensional case. In this case, the SSAP price should converge toward 
the theoretical arbitrage price when both the number of time steps d and the number of cells k 
converge towards infinity. Both European calls and European puts can be priced according to 
the original Black-Scholes formula. These prices are reported in the columns European Cbs 
and European Pgs of table (1). The American call can also be priced according to the same 
formula, since we assume the underlying asset pays no dividends. The price is reported in 
column American Cbs- F° r the American put, we computed the price using the finite-difference 
method presented in section (4.2). We call this method PDE, since it consists in solving a Partial 
Differential Equation. In dimension 1 , it is essentially equivalent to the Cox-Ross-Rubinstein 
binomial lattice method. We used 120 time steps for T (time to expiration) ranging from 1 to 4 
months, and 210 time steps for T - 7 months. The corresponding price is reported in column 
American Cpde- The SSAP prices where computed using M = 100, 000 samples, and k = 100 
buckets. The number of time steps was set to d = 10 in all the experiments. 

The observed differences between the SSAP prices and the reference prices are below 0.7%. 
The confidence interval values are below 1 % of the reference prices. American put prices given 
by the SSAP method are very accurate even when the difference with the European put prices 
are important (up to 30 cents). The computation time of a price using the SSAP method is about 
21 seconds, compared with less than one second with a classical integration method (PDE). In 
dimension 1 , classical finite difference of binomial lattice methods should be preferred to the 
SSAP method. 

8.3 Two underlying assets 

In this case the SSAP method only finds an approximation of the optimal price. However, 
numerical experiments show that the SSAP price always remains within a few cents of the 
optimal theoretical price. 
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Call option prices: x\ (0) = X2(0) = $40 








(Ji =20%, 02 : 


= 30%, r 


= 5% 




Parameters 


European 


American 


P 


1 




CpDE 




CpDE 


CsSAP 


H-j lilt, v 






35 


6.80 


6 79 0 010 

V -IS \J . V./ J V./ 


6.80 


6.80 


0.012 




1 

1 


4.0 


? 10 


9 11 0 00? 


? 10 


9 1 1 
z,. 1 1 


0 00? 
u.uu^ 






AS 


0 1 s 


0 18 0 00^ 


0 1 8 

U.lu 


0 1 8 

U. 1 0 


0 00^ 

U.UUJ) 






35 


8.90 


8 89 0 010 


8.90 


8.90 


0.010 


o <%. 

U /C 


A 


AO 


4 AS 


4 48 0 01 S 
'H-.'4-O U.U1J) 


4 48 


4 4Q 


0 01 S 
U.U1J) 






AS 


1 .00 


1 fifi 0 OOfi 
1.00 u.uuo 


1 .00 


1 .00 


0 008 
u.uuo 






35 


10.43 


10 42 0 020 

J V./ ■ 1 1mm V./ ■ V./ 1mm V./ 


10.43 


10.43 


0.020 




7 
/ 


4.0 


1 s 

U.l J 


1 S 0 016 

U.l J U.UIO 


fi 1 s 

O. LJ 


0. 1 o 


0 090 
u.u^u 






AS 


^ 1 0 


^ 08 0 004 
J.Uo U.UUM- 


^ 1 0 
J. 1U 


^ OQ 


0 010 
U.U1U 






35 


6.36 


6 35 0 004 


6.36 


6.36 


0.005 




1 
i 


40 


1 87 

1.0/ 


1 87 0 006 
1.0/ u.uuu 


1 87 

1.0/ 


1 87 
1 .o / 


0 008 
u.uuo 








0 1 7 
U.l/ 


0 17 0 004 

U.l/ U.ULH- 


0 1 7 
U.l/ 


0 1 7 
U.l/ 


0 004 
U.UUM- 






35 


8.09 


8.08 0.015 


8.09 


8.09 


0.016 


50% 


4 


40 


3.99 


3.98 0.010 


3.99 


3.99 


0.010 






45 


1.51 


1.50 0.006 


1.51 


1.51 


0.006 






35 


9.41 


9.40 0.015 


9.41 


9.41 


0.015 




7 


40 


5.48 


5.47 0.010 


5.48 


5.47 


0.016 






45 


2.77 


2.77 0.010 


2.77 


2.78 


0.012 






35 


5.61 


5.61 0.002 


5.61 


5.61 


0.006 




1 


40 


1.45 


1.45 0.008 


1.45 


1.46 


0.008 






45 


1.16 


1.16 0.002 


1.16 


1.16 


0.003 






35 


6.69 


6.68 0.003 


6.69 


6.69 


0.008 


100% 


4 


40 


3.09 


3.07 0.010 


3.09 


3.08 


0.012 






45 


1.24 


1.24 0.004 


1.24 


1.25 


0.010 






35 


7.62 


7.61 0.006 


7.62 


7.62 


0.012 




7 


40 


4.23 


4.21 0.020 


4.23 


4.22 


0.020 






45 


2.24 


2.22 0.012 


2.24 


2.23 


0.020 



Table 2: Prices for a call option with 2 underlying assets 
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Put option prices: x\ (0) = *2(0) = $40 


a y =20%,(72 = 30%,r = 5% 


Parameters 


European 


American 


P 


T 




Ppde 


Pssap 4stdev 


Ppde 


Pssap 4stdev 


0% 


1 


35 
40 
45 


0.00 
0.29 
3.34 


0 00 0 0000 
0.29 0.0030 
3.35 0.0060 


0.00 
0.40 
5.00 


0 00 0 0000 
0.40 0.0040 
5.00 0.0000 


4 


35 
40 
45 


0.03 
0.53 
2.62 


0 03 0 0020 
0.54 0.0080 
2.63 0.0200 


0.04 
0.75 
5.00 


0 04 0 0020 

v./ - v./ r v./ . v./ v./ j-. \j 

0.75 0.0080 
5.00 0.0000 


7 


35 
40 
45 


0.08 
0.66 
2.47 


0 08 0 0030 
0.67 0.0120 
2.46 0.0200 


0.10 
0.95 
5.00 


0 10 0 0040 

v./ - i v./ v./ - v./ v./ r v./ 

0.95 0.0100 
5.00 0.0000 


50% 


1 


35 
40 
45 


0.00 
0.50 
3.78 


0 00 0 0006 
0.50 0.0030 
3.78 0.0080 


0.00 
0.56 
5.00 


0 00 0 0006 
0.57 0.0050 
5.00 0.0000 


4 


35 
40 
45 


0.09 
0.92 
3.35 


0.09 0.0030 
0.91 0.0100 
3.35 0.0100 


0.10 
1.08 
5.00 


0.10 0.0030 
1.08 0.0100 
5.00 0.0000 


7 


35 
40 
45 


0.21 
1.14 
3.29 


0.21 0.0040 
1.13 0.0080 
3.30 0.0120 


0.24 
1.38 
5.00 


0.24 0.0040 
1.38 0.0120 
5.00 0.0000 


100% 


1 


35 
40 
45 


0.01 
0.83 
4.52 


0.01 0.0008 
0.83 0.0040 
4.52 0.0040 


0.01 
0.84 
5.00 


0.01 0.0006 
0.85 0.0060 
5.00 0.0000 


4 


35 
40 
45 


0.19 
1.51 
4.58 


0.19 0.0040 
1.51 0.0060 
4.59 0.0040 


0.19 
1.56 
5.02 


0.20 0.0040 
1.56 0.0120 
5.02 0.0120 


7 


35 
40 
45 


0.41 
1.87 
4.74 


0.41 0.0020 
1.86 0.0120 
4.73 0.0080 


0.42 
1.96 
5.20 


0.42 0.0060 
1.97 0.0100 
5.21 0.0120 



Table 3: Prices for a put option with 2 underlying assets 
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Call option prices: x\ (0) = x 2 (0) = X3 (0) 


= $40 








= 20%, a 2 = 30%, da = 50%, r = 5% 


Parameters 


European 


American 


P 


1 




CpDE 


CsSAP 




CpDE 


CsSAP 


H-Jttic- v 






35 


8.59 


8.58 


0.008 


8.59 


8.59 


0.010 




1 

1 


4.0 


^ 84 


^ 8^ 


0010 


% 8A 


^ 8A 


0 019 

U.U1Z 








O 8Q 

U.07 


O 8Q 

U.07 


O 006 


O 8Q 
U.oy 


O QO 
U.yU 


O 007 
u.uu / 






35 


12.55 


12.53 


0.020 


12.55 


12.55 


0.016 


o <%. 

U /C 


A 


AO 


7 87 
/ .0 / 


7 8^ 


O 01 A 


7 87 
/ .0 / 


7 87 
/ .0 / 


O 090 
U.UZU 






A*? 


A Ifi 


A 9^ 


O 01 A 


A Ifi 

H- .ZO 


A 97 
.Z / 


O 01 A 






35 


15.29 


15.27 


0.020 


15.29 


15.30 


0.030 




7 


4.0 


1 0 79 


1 0 70 


0 016 
u.u 1 u 


1 0 79 
iu. / z 


1 0 1% 


0 o^s 






A^ 


o.yo 


O.yJ 


O 090 

u.uzu 




6 Q8 


O 090 
u.uzu 






35 


7.78 


7.77 


0.010 


7.78 


7.78 


0.012 




1 
i 


40 


^ 1 8 


^ 17 


0 010 


^ 1 8 


^ 1 8 


0 019 

U.U1Z 






A 1 ^ 


O 89 
U.oZ 


O 89 


O OOA 
U.UUM- 


O 89 


O 8^ 


O OOA 
U.UUM- 






35 


10.97 


10.95 


0.010 


10.97 


10.96 


0.015 


50% 


4 


40 


6.69 


6.67 


0.016 


6.69 


6.69 


0.020 






45 


3.70 


3.69 


0.012 


3.70 


3.71 


0.025 






35 


13.23 


13.21 


0.016 


13.23 


13.24 


0.030 




7 


40 


9.11 


9.09 


0.020 


9.11 


9.12 


0.040 






45 


5.98 


5.98 


0.016 


5.98 


5.99 


0.020 






35 


6.53 


6.52 


0.006 


6.53 


6.54 


0.010 




1 


40 


2.38 


2.37 


0.010 


2.38 


2.38 


0.012 






45 


0.74 


0.74 


0.002 


0.74 


0.74 


0.003 






35 


8.51 


8.50 


0.008 


8.51 


8.53 


0.016 


100% 


4 


40 


4.92 


4.90 


0.012 


4.92 


4.93 


0.020 






45 


2.97 


2.96 


0.008 


2.97 


2.99 


0.020 






35 


10.04 


10.03 


0.010 


10.04 


10.07 


0.016 




7 


40 


6.64 


6.63 


0.016 


6.64 


6.67 


0.030 






45 


4.61 


4.60 


0.016 


4.61 


4.64 


0.040 



Table 4: Prices for a call option with 3 underlying assets 
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Put option prices: x\ (0) 


= * 2 (0) =x 3 (0) 


= $40 






0"! 


= 20%, a 2 = 30%, (73 = 50%, r = 5% 


Parameters 


European 


American 


P 


T 

1 




Ppde 


*SSAP ^siaev 


Ppde 


PsSAP 








35 


0.00 


0 00 0 000 

V./ - V./ V./ V./ - V./ V./ V./ 


0.00 


0.00 


0.000 




1 

1 


4.0 


0 1 % 

U. 1 J 


0 n 0 OfR 

U. LJ U.UUJ) 


0 1% 

U.x. J 


0 9^ 
u.z. j 


0 00^ 

U.UU J 






4.^ 

H-J 


9 96 


9 97 0 01 9 


S 00 


S 00 

J.UU 


0 000 






35 


0.01 


0 01 0 002 


0.01 


0.01 


0.001 


0 % 

u /c 


A 


40 


0 9S 




0 44 


0 4S 

U.T^ J 


0 OOfS 

U.UUU 






4.S 

H-J 


1 -J J 


1 Sfi 0 01 0 


S 00 

J.UU 


S 00 

J.UU 


0 000 

U.UUU 






35 


0.03 


0 03 0 004 


0.04 


0.04 


0.002 




7 


40 


0 ^1 

U. J 1 


0 ^9 0 01 0 


0 S7 
u. j i 


0 ss 

U.JO 


0 01 s 

U.U 1 J 






4.^ 

H-J 


1 41 


1 49 0 090 


S 00 

J.UU 


S 00 

J.UU 


0 000 
\J.\J\J\J 






35 


0.00 


0 00 0 000 

V./ - V./ V./ V./ - V./ V./ V./ 


0.00 


0.00 


0.000 




1 
i 


40 


o 

u. jo 


0 0 00^ 


0 4.8 


0 4.Q 


0 OOfi 

U.UUU 






4S 


^ 00 


^01 0010 


S 00 

J.UU 


S 00 

J.UU 


0 000 






35 


0.07 


0.08 0.006 


0.09 


0.09 


0.004 


50% 


4 


40 


0.72 


0.72 0.006 


0.93 


0.94 


0.010 






45 


2.65 


2.66 0.010 


5.00 


5.00 


0.000 






35 


0.17 


0.17 0.006 


0.20 


0.20 


0.005 




7 


40 


0.91 


0.91 0.012 


1.19 


1.21 


0.010 






45 


2.63 


2.65 0.012 


5.00 


5.00 


0.000 






35 


0.01 


0.01 0.001 


0.01 


0.01 


0.001 




1 


40 


0.84 


0.84 0.004 


0.84 


0.85 


0.006 






45 


4.18 


4.18 0.003 


5.00 


5.00 


0.000 






35 


0.19 


0.19 0.003 


0.19 


0.19 


0.004 


100% 


4 


40 


1.51 


1.51 0.006 


1.56 


1.57 


0.008 






45 


4.49 


4.49 0.005 


5.00 


5.00 


0.008 






35 


0.41 


0.41 0.004 


0.42 


0.42 


0.004 




7 


40 


1.87 


1.86 0.010 


1.96 


1.97 


0.012 






45 


4.70 


4.70 0.008 


5.20 


5.20 


0.012 



Table 5: Prices for a put option with 3 underlying assets 
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The European and American call and put option prices can be computed by the PDE method. 
This integration requires 120 time steps for T = 1 and T = 4 months and 210 time steps for 
T = 7 months. These results are reported in the columns European Cpde, Ppde and American 
Cpde, Ppde in tables (2) and (3). The SSAP method was run using M = 100, 000 samples, 
and k = 100 buckets. The number of time steps was set to d = 10 in all the experiments. 

The observed differences between the SSAP prices and the reference prices are below 1%. The 
confidence interval value is below 1% of the reference prices, except for very low prices where 
it remains under 1 cent. American put prices given by the SSAP method are very accurate. The 
computation time of a price using the SSAP method is about 25 seconds, compared with 23 
seconds for the classical integration method (PDE). In dimension 2, classical finite difference 
or binomial lattice methods are essentially equivalent to the SSAP method. 

8.4 Three underlying assets 

In this case again, the SSAP method only finds an approximation of the optimal price. However, 
numerical experiments show that the SSAP price always remains within a few cents of the 
optimal theoretical price. All parameters are identical to those of the 2D case, both for the PDE 
and the SSAP method. Results are presented in tables (4) and (5). 

The observed differences between the SSAP prices and the reference prices are below 1%. The 
confidence interval value is below 1% of the reference prices, except for very low prices where 
it remains under 1 cent. American put prices given by the SSAP method are very accurate. 
The computation time of a price using the SSAP method is about 32 seconds, compared with 
202 seconds for the classical integration method (PDE). In dimension 3, the SSAP method is 
as accurate and about 6 times faster than the classical integration method PDE. 

8.5 Ten underlying assets 

In the previous subsections 8.2, 8.3 and 8.4, we compared the efficiency and accuracy of 
the SSAP method with that of the classical integration method (PDE). In this subsection, we 
report results obtained with the SSAP method on American option pricing problems with 10 
underlying assets (tables (6) and (7)). Since no other method exists to compare to our results, 
and since the SSAP method only provides an approximation of the optimal price, we cannot 
guarantee the accuracy of the American premiums reported below. However, both the observed 
confidence intervals and the fact that the SSAP prices for the American calls without dividends 
equal the European prices lead us to believe that the SSAP method is reliable in general on 
10-dimensional American pricing problems. The parameters of the SSAP method are again 
M = 100, 000 and k = 100. 

The differences between European and American call prices are below 0.5 %. Since the payoff 
is the maximum of n underlying assets prices without dividends, these two prices should indeed 
be identical. Confidence intervals are around 1% in the worst case. The computation time 
using the SSAP method is about 82 seconds. 
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Table 6: Prices of a call option with 10 underlying assets 
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Put option prices:xi (0) = 
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Table 7: Prices of a put option with 10 underlying assets 
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Call confidence interval: Astdev 
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Put confidence interval: Astdev 
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Table 8: Efficiency of quadratic resampling (QR) 

8.6 Efficiency of extended quadratic resampling 

In this subsection we analyse the efficiency of the Quadratic Resampling method presented in 
section 7. The results shown in table (8) are the confidence intervals with (QR) and without 
(noQR) quadratic resampling. The gain value is the ratio of these two figures (noQR/ QR). We 
present in figure (1) the speed-up obtained through quadratic resampling. Since the accuracy 
of the Monte Carlo method is proportional to the square root of the number of samples (central 
limit theorem), the speedup obtained through Quadratic Resampling is the square of the gain 
(speed - up = gain 2 ). In these experiments, we used M = 10, 000 sample paths and k = 100 
buckets. 

Quadratic resampling is very efficient for a small number of underlying assets (up to 10). The 
speed-up ranges from 5 to 30, with an average of about 10. The speed-up decreases in higher 
dimensions. A speed-up of 2.5 is obtained with 40 underlying assets. 

8.7 Experimental time complexity 

We analyze in this subsection the computation time required by the SSAP method as a function 
of the number n of underlying assets. We compare these results when possible to those of 
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35 t 




Dimension 



Figure 1 : Speed-Up factor obtained through quadratic resampling 

classical integration methods (PDE). Observed computation times are presented in table (9). 
The number of Monte Carlo samples is M = 100, 000. The number of cells is k = 100 as 
before. 

Figure (2) shows that for up to 40 underlying assets, the computation time is linear in n. The 
dominating term is the computation of the payoff function which is linear in n (M x n x d). 
Figure (3) shows that for larger n the complexity is quadratic in n (M x n 2 x d), as expected. 

For n = 2 the SSAP is comparable to the classical integral method. For n = 3 the SSAP 
method is faster by a factor of 6. For n > 3 the SSAP method is the only one which can 
compute the price accurately. The integral (PDE) method is implemented using 60 time steps 
in the Cox-Rubinstein tree to obtain comparable precision. The exercise condition is also 
applied only 10 times during the life period of the option. 

8.8 Parallel implementation 

We implemented a parallel version of the SSAP method on a network of 4 DEC 3000 model 500 
ALPHA AXP workstations equipped with a high-b and width Gigas witch fiber optic interconnect 
(called a workstation farm). We observed a speedup linear in the number of workstations in 
the network: the parallel version is 4 times faster than the sequential implementation (table 
(10)). We anticipate that these figures would scale up with the number of workstations in 
the network. Indeed, the parallelization paradigm for Monte Carlo simulation is particularly 
simple. It consists in distributing on different processors the computation of different sample 
paths, and finally adding the results obtained by all processors. 

Let us assume that we have A independent computational units. The preliminary computations 
of the values a,-(f) , bij{t) , c,-(f) (see section 6) by Monte Carlo simulation can be done separately 
on each unit for Mj A Monte Carlo samples. Then they must be consolidated on one unit, called 
the master unit. The backward integration in time is then performed on the master unit using the 
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Table 9: Computation times as functions of the dimension 
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Figure 2: Linear behavior of computation time for 0 < n < 60 
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Table 10: Linear speedup of the parallel implementation 

consolidated values a,-(f) , btj(t) , Cj(t) . The memory requirement is multiplied by A as compared 
to a sequential implementation of the SSAP method. The gain in computation time is almost 
linear, i.e. the time complexity can be divided by A. Indeed, the results of the previous 
subsection have shown that the dominating terms are: (M x n x d) for the computation of 
the payoffs and (M x n 2 x d) for drawing the Monte Carlo sample paths. But these two 
operations are done on each computational unit separately. The computation time required for 
the backward integration can be neglected. (~ 1 second with the parameters of the previous 
subsection). The only overheads added by the parallelization are: 

• A communications of an amount of memory proportional to n 2 x d (typically 2 Mbytes). 

• Consolidation of the A sets of values. 

For n and d fixed, this overhead is constant and can be neglected in practice (~ 1 second for the 
parameters of the previous subsection). We present in table (10) observed figures for A = 4, 
and estimated figures for A = 32. 

9 Conclusion 

In this article, we described a systematic numerical technique for pricing arbitrarily complex 
American contingent claims, i.e. generalized option contracts with possibilities of early ex- 
ercise. Besides its obvious applications to trading and hedging in organized and Over The 
Counter (OTC) capital markets, American security pricing has many important applications 
in various areas of risk management such as assets and liabilities management and corporate 
investment decision making. Using this technique, we were able to compute the prices of 
complex American instruments in a few tens of seconds on a workstation, and within a few 
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seconds on a network of workstations. 

Our approach essentially relies on appropriate state aggregation techniques that circumvent 
the intractability of the computation of the early exercise boundary, combined with a classical 
Monte Carlo simulation for the computation of the conditional probabilities in the backwards 
pricing formula. We call this method Stratified State Aggregation along the Payoff function 
(SSAP). We have successfully implemented the SSAP method for problems with up to 400 
dimensions. To the best of our knowledge, no other method has ever been developed to date 
for pricing American contingent claims with many (more than 3 or 4) underlying assets. 

We feel that the method presented in this paper and the experimental results thus obtained 
make it possible to realistically envision the use of multidimensional stochastic models for 
practical real-world quantitative risk management problems. This capability of computing 
the joint influences of several tens of risk factors such as interest rates of various terms in 
different currencies, equity and commodities of various kinds, and any other relevant economic 
variables, may dramatically increase the competitive advantage of quantitative methods over 
more traditional analysis techniques. An application of particular interest is the pricing and 
hedging of complex long-dated commodity and index warrants offered on international OTC 
markets. We plan to backtest on actual market data the performance of the SSAP method as 
compared to more classical delta-hedging techniques currently used on capital markets. 
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